/**********************************************************************/
/*
   Title: analysis_FI_paper.do
   Authors: Robbie Dulin, Clotaire Boyer
   Created: 2020 11 24
   Description: Runs HH level analysis on SUSENAS 19
*/
/**********************************************************************/

cap log close
local prefix: display %tdCYND td(`c(current_date)')
log using "$log/`prefix'_HH_SUS", replace text

clear
set more off

set seed 11022022

global finance_data "$cleaned/finance"

** Globals
* global output                       "."
* global importable                   "."

qui cd "$output/ster_files/20210226_finance"

/*----------------------------------------------------*/
        /* Section: HH Regressions (Mar 19) */
/*----------------------------------------------------*/

  u "$finance_data/susenas_mar19_finance_treats_lasso.dta", clear

  * locals for regressions
  local 0    = "& num_agents_mar18 == 0"
  gen baseline_agent0 = num_agents_mar18 == 0
  local full = ""

local banking = "savings_account credit"
local pkh     = "pkh_receipt pkh_agent phk_other_location"

local cut_vars = "savings_account"

foreach var of varlist `banking' `pkh' {
  * get baseline control if it exists
  local mar18 = ""
  cap confirm variable `var'_mar18
  if _rc == 0 {
    local mar18 = "`var'_mar18"
  }
  foreach pmt in  all pmt30 above {

    foreach sample in 0 all {

      di "Var: `var' | PMT: `pmt' | Sample: `sample'"

      * DL regression
       { //quietly
        eststo `var'`pmt'`sample': pdslasso `var' treated (`mar18' mar16* mar17* mar18* *podes *udb i.finalstratum)  ///
          if `pmt' == 1 ``sample'', cluster(KABU) aset(i.finalstratum  `mar18')
          local effect = r(table)[1, 1]
          local se = _se[treated]
          local lower_bound = `effect' - 1.68 * `se'
          local upper_bound = `effect' + 1.68 * `se'

           summ `var' if e(sample) == 1 & treated == 0

          estadd scalar control_mean = r(mean)
          estadd scalar obs = e(N)
          estadd local stratum = "Yes"
          estadd local lasso = "Yes"
          // Calculate 90% confidence interval for treatment effect
          // which can be used in the paper (in-text)
      // Print confidence interval
      di "90% Confidence Interval for Treatment Effect: (`lower_bound', `upper_bound')"
      }

      *RI p-val
      randinference treated using "$importable/resample1000.dta", ///
        reg(reg `var' treated `e(xselected)' if `pmt' == 1 ``sample'', vce(cluster KABU)) ///
        reps($RI) mergevar(namakabupaten) mergetype(m:1)
      matrix define pvalues = (`r(RIpval)')
      matrix colnames pvalues = treated
      est restore `var'`pmt'`sample'
      estadd matrix pvalues
    }
  }
}



qui cd "$output/ster_files/20210226_finance"

* store estimates
  local prefix: display %tdCYND td(`c(current_date)')
  cap erase "`prefix'_banking_HH" // clear .ster if it already exists
  qui estimates dir
  foreach e in `r(names)' {
      quietly estimates restore `e'
      estimates title: `e'
      quietly estimates save "`prefix'_banking_HH", append
  }



  /*----------------------------------------------------*/
          /* Section: PKH GRAPH (Mar 19) */
  /*----------------------------------------------------*/

  u "$finance_data/susenas_mar19_finance_treats_lasso.dta", clear

  * locals for regressions
  local 0    = "& num_agents_mar18 == 0"
  gen baseline_agent0 = num_agents_mar18 == 0

  local pkh     = "pkh_receipt pkh_agent"

  foreach var of varlist `pkh' {
  * get baseline control if it exists
  local mar18 = ""
  cap confirm variable `var'_mar18
  if _rc == 0 {
    local mar18 = "`var'_mar18"
  }

    foreach sample in 0  {

      di "Var: `var' | PMT: `pmt' | Sample: `sample'"

      * DL regression
       { //quietly
        eststo `var'`pmt'`sample': pdslasso `var' treated (`mar18' mar16* mar17* mar18* *podes *udb i.finalstratum credit_mar18)  ///
          if pmt30 == 1 & num_agents_mar18 == 0, cluster(KABU) aset(i.finalstratum  `mar18')
          local treatment_effect_`var' = r(table)[1, 1]
          local se_`var' = _se[treated]
          summ `var' if e(sample) == 1 & treated == 0
          estadd scalar control_mean = r(mean)
          estadd scalar obs = e(N)
          estadd local stratum = "Yes"
          estadd local lasso = "Yes"

      }

    }
  }


// We want to create a plot for the PKH results

  * Compute the mean for treatment group
  summarize pkh_receipt if treated == 1 & pmt30 == 1 & num_agents_mar18 == 0, meanonly
  local treatment_mean_receipt = r(mean)
  summarize pkh_agent if treated == 1 & pmt30 == 1 & num_agents_mar18 == 0, meanonly
  local treatment_mean_agent = r(mean)

  * Compute the control value as the mean minus treatment effect
  local control_value_receipt = `treatment_mean_receipt' - `treatment_effect_pkh_receipt'
  local control_value_agent = `treatment_mean_agent' - `treatment_effect_pkh_agent'

  * Compute the confidence interval for the difference in means
  local lower_ci_receipt = `control_value_receipt' - 1.645 * `se_pkh_receipt'
  local upper_ci_receipt = `control_value_receipt' + 1.645 * `se_pkh_receipt'

  local lower_ci_agent = `control_value_agent' - 1.645 * `se_pkh_agent'
  local upper_ci_agent = `control_value_agent' + 1.645 * `se_pkh_agent'



  * Create a matrix to store values
  matrix pkh_values = J(4, 5, .)

  * Extract numeric values from local macros and assign to matrix elements
  matrix pkh_values[1, 1] = real("`treatment_mean_receipt'")
  matrix pkh_values[2, 1] = real("`treatment_mean_agent'")
  matrix pkh_values[3, 1] = real("`control_value_receipt'")
  matrix pkh_values[4, 1] = real("`control_value_agent'")

  * Compute confidence intervals and store in matrix
  matrix pkh_values[3, 2] = real("`lower_ci_receipt'")
  matrix pkh_values[3, 3] = real("`upper_ci_receipt'")
  matrix pkh_values[4, 2] = real("`lower_ci_agent'")
  matrix pkh_values[4, 3] = real("`upper_ci_agent'")

  * Compute confidence intervals and store in matrix
  matrix pkh_values[1, 4] = 1
  matrix pkh_values[2, 4] = 2
  matrix pkh_values[3, 4] = 1
  matrix pkh_values[4, 4] = 2

  * Compute confidence intervals and store in matrix
  matrix pkh_values[1, 5] = 1
  matrix pkh_values[2, 5] = 1
  matrix pkh_values[3, 5] = 2
  matrix pkh_values[4, 5] = 2

  * Save matrix as dataset
  clear
  svmat pkh_values
browse
rename pkh_values1 value
rename pkh_values2 ci_moins
rename pkh_values3 ci_plus
rename pkh_values4 outcome
rename pkh_values5 group

* Create a variable to group Treatment 1 and Treatment 2
gen outcomegroup = outcome
replace outcomegroup = outcome + 2.4 if outcome > 1
replace outcomegroup = outcomegroup + 1.2 if group == 2

* Bar plot with gaps
graph twoway (bar value outcomegroup if group==1,color(maroon) fintensity(inten60)) ///
 (bar value outcomegroup if group==2,color(navy) fintensity(inten60)) ///
 (rcap ci_plus ci_moins outcomegroup), ///
  legend(order (1 "Voucher Mean" 2 "Counterfactual Mean" 3 "CI - 90%")) ///
 xlabel(1.6 "PKH Receipt" 4.95 "PKH Delivered by Bank Agent", noticks) ///
 ytitle("Proportion of Beneficiaries") ///
 yscale(range(0 1)) ylabel(0 0.2 0.4 0.6 0.8 1) xtitle("")

graph export "$output/figures/pkh.png", replace


cap log close

// DONE
